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Filename: SPIRALI.FOR 


c *********************************************************************** 

C UPDATED INTERIM VERSION OF SPIRALI AS OF 3/9/95 THAT INCLUDES 
C LOCAL PRESSURE JUMPS AT GROOVE-RIDGE INTERFACES. ALTHOUGH THIS 
C COOE DOES PROVIDE INPUT FOR CARRYOVER EFFECTS IT IS INCOMPLETE IN THAT 
C I t ^DOES^WO T ^FULL Y ^IM P LEHEWT^TH H M^AS^OF^YET.^ aa j^^ a ^ a ^ x 11JLJL- L 1J L 1Jl J lJUJU 

C MAIN PROGRAM FOR COMPUTER COOE SPIRALJ WHICH IS EXTENDS SPIRALI TO 
C INCLUDE LOCAL PRESSURE JUMPS AT GROOVE-RIDGE DISCONTINUITIES 
C FL /Gt1024 SPIRALI.FOR 

PARAMETER <NDZ=201,NDREG=21) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

CHARACTER T I TLE*64 , PNAME*60, F NAME (3 )*60 
DIMENSION NRSUB(NDREG),ELFR(NDREG),RH20(2) 

DIMENSION ALP I (NDREG) , BET I ( NDREG ) ,DELT ( NDREG) , DELT I ( NDREG ) 

DIMENSION NSG(NDREG) , ENGP(NDREG) , ZETG(NDREG) 

DIMENSION RG(NDZ, NDREG), ZET( NDREG ),ZTG(NDZ, NDREG) 

DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ,NDREG),P(NDZ, NDREG) 

DIMENSION TAU(NDZ,NDREG),CK(4,4),CB(4,4),AM<4,4) 

C THIS COMMON BLOCK PASSES DATA TO USER DEFINABLE FUNCTION FLMSHP 
COMMON/BFSHP/HTAP , HBRL 

C THIS COMMON BLOCK PASSES TURBULENCE COEFFICIENTS TO 
C USER DEFINABLE FUNCTIONS FA AND FB 
COMMON/BFAFB/EMA, ENA, EMB.ENB 
DATA RH20/9. 357260-5,1. 03/ 

C INITIALIZE NAMELIST VARIABLES 

DATA TITLE, I FACE, ISIUN, IGROT.NOI , IFLOW/' ',5*0/, 
+R0,C,EL,RPM,RPM0,RPMD/6*0.D0/, 

+PLEG, PRIG, VISC, DENS/4*0 . DO/ , 

+FZD, IH0ME,NITH,T0LH/0. DO, 0,10,1 .D-4/ 

+NITV,TOLV,DUT/30, 1 .D-5, 1 .0-6/ , 

+NREG,NRSUB(1),ELFR(1)/l,20,1 .DO/ , * n , 

+ALPI ,BETI ,DELT,ZET/NDREG*O.DO,NDREG*O.DO,NDREG*O.DO,NDREG*O.DO/ , 

+NSG, ZETG/NDREG*0 , NDREG*0 .DO/ 

NAMELIST/INPUTS/TITLE, IFACE, ISIUN, IGROT.NOI, IFLOW, 

+R0, C, EL, RPM,RPMO,RPMD,PLEG, PRIG, VISC, DENS, EMA, ENA, EMB.ENB, 

♦HTAP , HBRL , FZD , I HOME , N I TH , TOLH , 

♦NITV.TOLV.DUT, 

+NREG,NRSUB,ELFR,ZET,ALPI,BETI,DELT,NSG,ZETG „,cn 

C USE DATA STAEMENT BELOW TO HARD CODE DEFAULT FILENAMES. BLANK VALUES USED 
C HERE CAUSES MICROSOFT COMPILER TO TAKE NAMES FROM COMMAND LINE OR ISSUE 
C PROMPTS FOR THEIR INPUT AT RUN TIME. 

DATA FNAME/3*' '/ 

C UNIT 1 IS INPUT FILE, 2 OUTPUT FILE, 3 PLO T FI LE 
OPEN < 1 , F I LE= FNAME < 1 ) , ST ATUS= ' OLD ' , ERR=9999 ) 
0PEN(2,FILE=FNAME(2),ERR=9999) 

0PEN(3,FILE=FNAME(3),ERR=9999) 

INQUIRE(3,NAME=PNAME) 

PI=4.D0*ATAN(1 .DO) 

I CASE-0 

C INITIALIZE DATA ON COMMON BLOCKS 
HTAP=O.DO 
HBRL=O.DO 
EMA=-0.25D0 
EMB=EMA 
ENA=.0791D0 
ENB=ENA 

1 READ(1, INPUTS, END*999) 

ENGP b NSG/4 . DO/PI 
C CLEAN UP FLAGS 

I F ( I GROT . NE . 1 ) I GROT=0 
IF(IFACE.NE.1)IFACE=0 

I F( ISIUN. NE.1) ISIUN >0 ...... 

C SET CASE NUMBER AND WRITE CASE NUMBER AND TITLE TO OUTPUT FILE 
1CASE=ICASE+1 

IF(ICASE.GT.1)WRITE<2,*)' ' 


1_ 04/14/1995 13:06 Filename: SPIRALI.F OR ES9S — Z 

WRITE(2,19)ICASE, TITLE 
19 FORMATS ( CASE', 13,' ) ',A64/) 

C PRINT OUT NAMELIST . 

CALL INLIST(2,TITLE, IFACE, ISIUN, IGROT.NOI, IFLOW, 

+R0, C, EL, RPM,RPMO,RPMD,PLEG, PRIG, VISC, DENS, 

+EMA,ENA,EMB , ENB, HTAP , HBRL , 

♦FZD, IHOME.NITH, TOLH, 

+NITV,TOLV,DUT, 

+NREG,NRSUB,ELFR,ZET,ALPI,8ETI,DELT,NSG,ZETG) 

CALL OUTSCR(' STARTING SOLUTION FOR CASE NUMBER', I CASE) 

C CHECK UP FRONT FOR ERRORS 

IF(RO.LE.O.DO.OR.EL.LE.O.DO.OR.C.LE.O.DO.OR.VISC.LE.O.DO 
+.0R.RPM.LT.-1 .D-8)IER=8 
IF(NREG.GT.NDREG)IER=10 
ELSUM=O.DO 
DO 81 K*1,NREG 
ELSUM*ELSUM+ELFR(K) 

IF<NR$UB(K).GE.NDZ)IER*9 

81 CONTINUE 

I F( ABS( ELSUM- 1 .DO) . GT . 1 . 01 ) I ER*1 1 
IF(IER.NE.O)GO TO 80 
C GENERATE R AND Z GRIDS 
ELT»EL/2.D0/R0 

CALL RZGRIDt IFACE, ELT,NREG,NRSUB,ELFR,RG,ZTG) 

C CHECK ON HOMING IN ON AXIAL LOAD FOR FACE SEAL 
ITH*0 

IF(IHGME.NE.2)CNEW*C 

IF(IFACE.EQ.1 .AND.(IH0ME.EQ.1.0R.IH0ME.EQ.2) 

+.AND.FZD.GT.1 .D-6)ITH«1 

C LABEL BELOW IS TOP OF HOMING LOOP, USED WHEN ITH>0 
88 C1=CNEW 

IF(ITH.GT.O)CALL OUTSCR(' LOAD ITERATION NO.'.ITH) 

C GET rIfERENCE GAGE PRESSURE AND CHECK DIRECTION OF POISEUILLE FLOW 
PO=PLEG-PRIG 
IDIR=1 

IF(PO.LT.O.)IDIR=-1 
PO=ABS(PO) 

C CALCULATE DIMENSIONLESS PRIMARY FILM THICKNESS H, AND TOTAL NO. PTS, MTOT 
MTOT»0 

DO 7 K"1,NREG 
NRT*NRSUB(K)+1 
MTOT«MTOT+NRT 
DELTI (K)=DELT (K)/C1 
ENGP(K)*NSG(K)/4. DO/PI 
DO 7 J-1.NRT 

C ADD SHAPE AND DIVIDE BY C 
X=ZTG(J,K)/(2.D0*ELT) 

IF(IFACE.EQ.1)X b X+(ELT-1.D0)/(2.D0*ELT) 

7 H(J,K)*1 .DO+DELTI(K)*ALPI (K)+FLMSHP(X)/C1 

C USE SMALL NUMBER IN PLACE OF 0 DENSITY 
DNS=MAX< 1 .D-8*RH20( ISIUN+1 ) ,DENS) 

C CONVERT ANGULAR VELOCITIES TO RAD/SEC 
DRC=PI/30.D0 
OM=RPM*DRC 
OMO=RPMO*DRC 
OMD=RPMD*DRC 

C CALCULATE VELOCITY AND REYNOLDS NUMBER FOR LAMINAR FLOW 
VPL=C1*C1/(12.D0*VISC*EL) 

VL*VPL*PO 

REL=2 . D0*C1 *VL*DNS/VI SC 
REA=REL 
RFBR=24 DO 

C GET CHARACTERISTIC AXIAL REYNOLDS NUMBER FOR TURBULENT FLOW 
IF(REA.GT.1000.D0)THEN 
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CALL RECAL(REL,REA,FBAR,NITV,DUT , IER) 

IF(IER.NE.O)GO TO 99 
RFBR=REA*FBAR 
END IF 

VPT =VPL*24 . DO/RFBR 

C COMPARE AXIAL AND CIRCUMFERENTIAL REYNOLDS NUMBERS TO SELECT RE AND PO 
RE0=C1 *RO*OM*DNS/V I SC 
I F ( REA. GE . REO)THEN 
RE=REA 
PIN=1.D0 
ELSE 
TMP=PO 

PO=RO*OM/2./VPT 
PIN=TMP/PO 
RE=REO 
END IF 

C MAKE SOME DECISIONS ON HANDLING AXIAL INERTIA 
N0I1-0 

IF(NOI .GE.0.AND.(C1/EL*REA/RFBR.LT . .01 .OR.NOI .GT.O))NOI1*1 
IF(IFLOW*IDIR.EQ.-1 .AND.NOI1 .EQ.0)THEN 
IDIR--IDIR 
PIN=-PIN 
END IF 

REC=RE*2.D0*C1/R0*DENS/DNS 

IF(NOI.EQ.2)REC«O.DO 

C COMPUTE CHARACTERISTIC VELOCITY AND DIMENSIONLESS PARAMETERS 
V0=RE*VI SC/ < 2 . D0*DNS*C1 ) 

P1R=VISC*V0*R0/(4.D0*C1*C1*P0) 

OMT=OM*RO/VO 
OMOT =OM0*R0/V0 
OMDT =OMD*R0/V0 

C CALCULATE DIMENSIONLESS LOADING 
FZND=FZD/P0/R0**2 
C PERFORM SEAL COMPUTATIONS 

CALL TSEAL( TOLV , N I TV, NOI 1 , I FACE , I D I R , I GROT , NREG , NRSUB , 
+RE,REC,P1R,PIN,0MT,0M0T,0MDT,DUT,ZET,ALPI ,BETI ,DELTI ,ENGP,ZETG, 
+RG,ZTG,H,U,V,P,TAU,CK,CB,AM,FLO,TOR,W, I AMASS, ITER, IER) 

PADD=0 . DO 

IFCPIN.LT. O.DO)PADD=-PIN 
JEND=NRSUB(NREG)+1 

V=W+PI*PADD*(RG(JEND,NREG)**2-RG(1,1)**2) 

RE1=RE/DNS*DENS 

REA1 =RE1*H( 1,1 )*ABS(V< 1,1)) 

REA2=RE1*H( JEND,NREG)*ABS(V(JEND,NREG)) 

REOl=REl*H(1 , 1 )*ABS(U( 1 , 1 )-RG( 1 , 1 )*OMT) 

RE02=RE1*H(JEND,NREG)*ABS(U(JEND,NREG)-RG(JEND,NREG)*0MT) 

C CHECK TO SEE IF RIGHT INLET BOUNDARY WAS USED 

IFCIER.EQ.O.AND.NOI1 .EQ.O.AND.IFLOW*FLO.LT.O. )IER=7 
C IF HOMING ON LOAD FOR FACE SEAL, CHECK FOR CONVERGENCE OR TROUBLE 
IFCITH.GT.O.AND.IER.EQ.O.JTHEN 
I FCABSCW- FZND )/FZND . LT . TOLH )THEN 
CONTINUE 

ELSE IF(ITH.EQ.NITH)THEN 
IER=5 

ELSE IF(CK(1,1).LT.1.D*20)THEN 

IER=6 

ELSE 

CNEW*( 1 .DO- ( FZND-W)/CK(1 , 1 ) )*C1 
I FCCNEU/C1 . LT . 1 .D-8)THEN 
IER=6 
ELSE 

ITH=ITH+1 
GO TO 88 
END IF 
END IF 
END IF 
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C WRITE OUTPUT 
N0I2=N0I 1 

I FCNOI 1 .EQ.1 .AND.REC.LT. 1 .D-4)N0I 2=2 
80 CALL DIMOUT(2,NOI2,IFACE,ISIUN, IER, ITER, 

+W FLO TOR CIC CB AM 

+r6,el!ci ,visc,dens'pleg,prig,vo,po, iamass,rpm,rpmo,rpmd, 
+REA1,REA2,RE01,RE02) 

c if non-nul"pl6? ) f?le°dump no. points, film, velocity and pressure data 

IFCPNAME.EQ. 'NUL' .OR.PNAME.EQ. 'nul' )GO TO 99 
WRITEC3 ,*)MTOT 
ROM2=OM*RO/2.DO 
DO 20 K=1,NREG 
NRT=NRSUB(K)+1 

20 WRITE(3,21)(ZTG(I,K),H(I ,K)*C1-ALPI(K)*DELT(K), 
+U(I,K)*VO,V(I,K)*VO,(P<I,K)+PADD)*PO, I=1,NRT) 

C 21 +U??^K)*VO|v(l)K)*v6,(P<i,K)+PADD)*PO,RG(I,K)*ROM2, 1*1 ,NRT) 

C 21 FORMAT (OP, FI 0.4, IP, 5E13. 5) 

99 IF(IER.NE.0)CALL EMSG(IER) 

WRITEC2,*)' ' 

GO TO 1 
999 CLOSE(I) 

CLOSEC2) 

CLOSEC3) 

9999 STOP 
END 
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SUBROUTINE EMSG(IER) 

CHARACTERV8 MSG(11) 

C SENDS ERROR MESSAGES TO STD. OUTPUT. 

C CALLED BY MAIN PROGRAM 
DATA MSG/ 

♦'INITIAL VELOCITY COMPUTATION DIVERGED', 

♦'PRIMARY FLOW COMPUTATION DIVERGED', 

♦'MATRIX INVERSION ERROR ENCOUNTERED IN SECOND ORDER SOLUTION', 
♦'SPIRAL GROOVE LOCAL FLOW COMPUTATION DIVERGED', 

♦'FACE SEAL AXIAL LOAD ITERATION DIVERGED', 

♦'NEGATIVE STIFFNESS OR FILM THICKNESS IN AXIAL LOAD ITERATION', 
♦'WRONG INLET BOUNDARY WAS USED WITH TRANVERSE INERTIA INLCUDED' , 
♦'ILLEGAL LENGTH, CLEARANCE, VISCOSITY, PRESSURE OR SPEED ENCOUNTER 
♦ED', 

♦'MAXIMUM NUMBER OF ALLOWABLE GRID POINTS EXCEEDED', 

♦'MAXIMUM NUMBER OF ALLOWABLE REGIONS EXCEEDED', 

♦'SUM OF LENGTH FRACTIONS ARE NOT EQUAL TO 1' 

+/ 

WRITER, *>MSG(IER) 

RETURN 

END 
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SUBROUTINE OUTSCR(MSG,NUM) 

C SENDS STATUS MESSAGES TO THE STANDARD OUTPUT UNIT 
C CALLED BY MAIN, TSEAL 

CHARACTER* ( * )MSG , CNUM*6, MSG 1 *78 
WRITE<CNUM,'<I6)')NUM 

DO 5 11=1,6 
1=11 

IF(CNUM(I:I).GT.' ') GO TO 6 

5 CONTINUE 

6 CONTINUE 
L=LEN(MSG) 

MSG1=MSG 

C CONCATINATE NON 0 NUMBER TO STRING 
IF(NUM.EQ.O)THEN 
LT=L 
ELSE 
LT=L+8-I 

MSGKL+1 :LT) =/ '//CNUM<I:6) 

END IF 

WRITE(*,*)MSG1(1 :LT) 

RETURN 

END 
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SUBROUTINE DIMOUT(IFILE,NOI,IFACE,ISIUN,IER,ITER, 
+W,QIN,TOR,CK,CB,AM, 

+DDR , DDL, DDC , DDMU , PENS , DDPL , DDPR , VO , PO , I AMASS , RPM , RPMO , RPMD , 
+REA1,REA2,RE01,RE02) 

C SENDS OUTPUT TO UNIT I FILE 
C CALLED BY MAIN 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

CHARACTER STRK(4,2)*7,STRK0(4,2)*7.STRB(4.2)*7, 
+STRA(4,2)*7,KUNIT(4,2.2)*10,BUNIT(4,2,2)*j0,AUNIT(4,2,2)*10 
CHARACTER IN(2)*6.LB(2)*6, I3S<2)*13 J .PSI(2)*7. IHP(2)*8, 
+XIN(2)*12,YIN(2)*12.ZIN(2)*12,ILBC(2)*10,PSI4(2)*14,PSIS(2)*12 
CHARACTER NOISTR(3)*30 

DIMENSION CK(4,4),CB(4,4),AM(4,4),SC0N(4),XC0N(4) 

DATA NOISTR/' ',', TRANSVERSE INERTIA NEGLECTED', 

+', INERTIA NEGLECTED'/ 

DATA STRK/'Kx','Ky' 'Kphi','Kpsi' 'Kz','ICphi','Kpsi',' V ## 

DATA STRKO/'KOx' , 'KOy' , 'KOphi ' , 'KOpsi ' , 'KOz' , 'KOphi ' , 'KOpsi ' , ' '/ 
DATA STRB/'Bx' ,'By' , 'Bphi' , 'Bpsi ' , 'Bz' , 'Bphi ' , 'Bpsi ', ' '/ 

DATA STRA/'Ax','Ay','Aphi','Apsi','Az','Aphi','Apsi',' '/ 

DATA KUNIT/'LB', 'LB', 'IN-LB', 'IN-LB', 'LB', 'IN-LB', 'IN-LB',' ', 
+'N','N','N-m','N-m','N','N-m','N-m',' '/ 

DATA BUN IT/ 'LB- SEC' , 'LB-SEC' , ' IN-LB-SEC' , 'IN-LB-SEC' , 'LB-SEC' , 

+ MN-LB-SEC', 'IN-LB-SEC',' ', 

+'N-SEC' , 'N-SEC' , 'N-m-SEC' , 'N-m-SEC' , 'N-SEC' , 

+ DATA AUN I T / ' LB - slc2 ' , ' LB - SEC2 ' , ' IN-LB-SEC2' , ' IN-LB-SEC2' , 

+' LB-SEC2' , ' IN-LB-SEC2' , ' IN-LB-SEC2' , ' ' , 

+'N-SEC2' 'N-SEC2' , 'N-m-SEC2' , 'N-m-SEC2' , 'N-SEC2' , 
+'N-m-SEC2','N-m-SEC2',' '/ 

DATA XIN,YIN,ZIN/' x (IN)',' x (m)',' y (IN)',' y (m)', 

+' z (IN)',' z (m)'/ 

DATA IN, LB/' (IN)',' (m)',' (LB)',' (N)'/ 

DATA ILBC/' (IN-LB) ',' (N-m),'/ 

DATA I3S,PSI/' (IN**3/SEC)' , ' (m**3/SEC)' , ' (PSI)',' (Pa)'/ 

DATA PSIS/' (PSI-SEC),',' (Pa-SEC),'/ 

DATA PSI4/' (LB-SEC/IN4)',' (Kg/m3)'/ 

DATA I HP/' (HP)',' (WATT)'/ 

K=1 

IF( ISIUN.EQ.1 )K=2 
IFdFACE.EQ.DGO TO 1000 
IF(IER.EQ.O)THEN 
FCON=PO*DDR*DDR 
QCON=VO*DDC*DDR 
TCON=DDR/VO 

HP=T0R*FC0N*DDC*RPM*1 .586662957D-5 
I F ( X . EQ . 2 ) HP=HP*6600 
SC0N(1 )=FCON 
SC0N(2)=FC0N 
SC0N(3)=F CON*DD R 
SC0N(4)=FC0N*DDR 
XC0N(1 )=DDC 
XC0N(2)=DDC 
XC0N(3) S DDC/DDR 
XC0N(4)=DDC/DDR 
END IF 
NFP=4 

WRITE(IFILE,60)NOISTR(NOI+1 ) 

60 FORMAT (/' CYLINDRICAL SEAL',A30/) 

WR I TE ( I F I LE , 61 )DDL,2.D0*DDR,DDC, IN(K) „ , 

61 FORMAT ( ' LENGTH, DIAMETER, CLEARANCE =',1P,E12.4,',',E12.4,',', 
+E12.4,A6/) 

WRITE(IFILE,62)RPM, RPMO, RPMD 

62 FORMAT ( ' ROTOR, SWIRL AND DIST. SPEEDS *',1P,E12.4, 

+' , ' ,E12.4, ' ,' ,E12.4, ' (RPM)'/) 

WRITE(IFILE,63)DDPL,DDPR,PSI(K) „ ^ , 

63 FORMAT ( ' PRESSURE AT START, END AXIAL BOUNDARIES *' ,1P,E12.4, ' , ' , 



64 

40 

48 

50 

45 

100 

101 


+E12.4,A7/) 

WRITE(IFILE,64)DDMU,PSIS(K) 4 DENS.PSI4(K) 

FORMAT ( ' VISCOSITY =' , 1P,E12.4,A12, 

+' DENSITY *' ,E12-4,A14/) 

WRITE(IFILE,40)IER,ITER 

IF(IER.NE.O)RETURN , # 

FORMAT ( ' ERROR CODE =',I3 ', ITERATIONS IN PRIMARY FLOW =',I3) 
WRITE(IFILE,48)QIN*QCON,l4s(K),TOR*FCON*DDC,ILBC(K), 

+HP, IHP(K),REA1 REOl ,RE02 
FORMAT (/' FLOW *' 1P,E12.4,A13// 

+' TORQUE *',1P,El5.4,A10, ' FILM POWER LOSS =' ,E12.4,A8// 

+' AXIAL REYNOLDS NUMBER =' ,1P,E12.4,/ 

+' CIRC. REYNOLDS NUMBERS FOR ROTOR AT SEAL ENDS =' , 1P,E12.4, ' , ' , 


FORCE UNIT 


+E12.4) 

WRITE( I FILE,50) ^ 

FORMAT(/' DYNAMIC COEFFICIENTS ( FORCE UNIT / DISP. UNIT 

+)') 

ASSIGN 45 TO KF 
WRITE(IFILE,KF)XIN(K),YIN(K) 

FORMAT (/' DISP. ' 2A12, 

+' phi (RAD) ',' psi (RAD) ',' FORCE UNIT ') 

ASSIGN 47 TO KF 
DO 100 1*1, NFP 

I F( IAMASS.EQ.1 )THEN 
WR I TE( I F I LE , KF )STRK( I , I FACE+1 ) , 

+ (CK(I , J)*SCON(I )/XCON( J), J=1 ,NFP),KUNIT(I , IFACE+1 ,K) 

ELSE 

WR I TE( I F I LE , KF )STRK( I , I FACE+1 ) , 

+ (AM(I , J)*SCON(I )/XCON( J), J=1 ,NFP),KUNIT(I , IFACE+1 ,K) 


CONTINUE 
DO 101 1*1, NFP 

WRITE( I FILE,KF)STRB( I , I FACE+1 ) , (CB( I , 
+J*1. NFP), BUNIT(I, IFACE+1, K) 

DO 102 1*1, NFP 

IF(IAMASS.EQ.1)THEN 
WR I TE(I F I LE , KF )STRA(I , I FACE+1 ) , 

+ (AM( I , J )*SCON ( I )*TCON**2/XCON (J), 

ELSE 

WR I TE ( I F I LE , KF ) STRKO (1,1 FACE+1 ) , 

+ (CK( I , J )*SCON( I )/XCON( J ) , J*1 ,NFP) 

END IF 


J )*TCON*SCON ( I )/XCON (J), 

J*1,NFP),AUNIT(I, IFACE+1,K) 
,KUNIT(I, IFACE+1, K) 


102 CONTINUE 

47 FORMAT(1X,A7,1P,4E12.4,3X,A10) 

RETURN 

1000 IF(IER.EQ.O)THEN 
FCON*PO*DDR*DDR 
QCON*VO*DDC*DDR 
TCON*DDR/VO 

HP*TOR*FCON*DDC*RPM*1 .586662957D-5 
IF(K.EQ.2)HP*HP*6600 
SCON(1)*FCON 
SCON(2)*FCON*DDR 
SCON(3)*FCON*DDR 
XCON(1)*DDC 
XCON ( 2 )*DDC/DDR 
XCON(3)=DDC/DDR 
END IF 
NFP*3 

WRITE(IFILE,860)NOISTR(NOI+1) 

860 FORMAT (/' FACE SEAL',A30/) 
WRITE(IFILE,861)2.D0*(DDR-DDL),2.D0*DDR,DDC,IN(K) 

861 FORMAT (' ID, 00, NOMINAL FILM THICKNESS *', 
+1P,E12.4,',',E12.4,',' E12.4,A6/) 

WRITE(IFILE,62)RPM, RPMO, RPMD 
WRITE(IFILE,863)DDPL,DDPR,PSI(K) 
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863 FORMAT ( ' INSIDE, OUTSIDE PRESSURE =' , IP, E12.4, ', ' , 

+E12.4.A7/) 

URITEd FILE, 64)DDMU,PSIS(K>, DENS, PSI4(K) 

WRITE(IFILE,40)IER,ITER 

IF(IER.NE.O)RETURN 

URITEdFILE, 144)U*FCON,LB(K) „ . 

144 F0RMAT(/' AXIAL LOAD TO BALANCE FACE SEAL =' ,1P,E12.4,A6) 

URITEd FILE, 848)QIN*QCON,I3S(K),TOR*FCON*DDC,ILBC(K), 

+HP,IHP(K),REA1, REA2 , RE01 , RE02 
848 FORMAT(/' FLOW = ' 1P,E12.4,A13// 

+' TORQUE =',1P,E12.4,A10,' FILM POWER LOSS =' .E12.4.A8// 

+' RADIAL REYNOLDS NUMBER AT ID, OD =',1P,E12.4 ,',',E12.< 

+/' CIRC. REYNOLDS NUMBERS FOR ROTOR AT ID, OD =',1P,E12.4,',', 
+E12.4) 

URITEd FILE, 50) 

ASSIGN 145 TO KF 
URITEd FI LE,KF)ZIN(K) 

145 FORMATS' DISP. ' A12, 

+' phi (RAD) ',' psi (RAD) ',' FORCE UNIT ') 

ASSIGN 147 TO KF 
DO 500 1=1, NFP 

IF(IAMASS.EQ.1)THEN 

URITE(IFILE,KF)STRK(I,IFACE+1), 

+ (CK(I , J)*SCON(I )/XC0N(J), J=1 ,NFP),KUNIT(I, IFACE+1 ,K) 

ELSE 

UR I TE (I F I LE , KF )STRK(I , I FACE+1 ) , 

+ (AM(I,J)*SCON(I)/XCON(J),J=1, NFP), KUNITd, IFACE+1, K) 

END IF 

500 CONTINUE 

DO 501 1=1, NFP 

501 URITEd FILE, KF)STRB(I, I FACE+1 ),(CB(I , J)*TCON*SCON( I )/XCON(J), 

+J=1, NFP), BUNIT( I, IFACE+1, K) 

DO 502 1=1, NFP 

I F ( I AMASS . EQ . 1 )THEN 

URITEdFILE, KF)STRA(I, IFACE+1), 

+ (AM(I,J)*SCON(I)*TCON**2/XCON(J),J=1,NFP),AUNIT(I,IFACE+1,K) 

ELSE 

UR I TE(I F I LE , KF )STRK0( I , I FACE+1 ) , 

+ (CK(l,J)*SCON(I)/XCON(J),J=1, NFP), KUNITd, IFACE+1, K) 

END IF 

502 CONTINUE 

147 FORMAT(1X,A7, 1P,3E12.4,3X,A10) 

RETURN 

END 


SUBROUTINE INLIST( IFILE, TITLE, I FACE, ISIUN, IGROT,NOI , I FLOU, 
+RO,C,EL,RPM,RPMO,RPMD,PLEG,PRIG,VISC,DENS, 

+EMA , ENA, EMB, ENB , HTAP , HBRL , 

+FZD, IHOME,NITH,TOLH, 

+NITV,TOLV,DUT , 

+NREG,NRSUB,ELFR,ZET,ALPI,BETI ,DELT,NSG,ZETG) 

C THIS ROUTINE PRINTS OUT THE NAMELIST IN A LEGIBLE MANNER 
C CALLED BY MAIN PROGRAM 
C WRITES TO UNIT NO. IFILE 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

PARAMETER (NDZ=201,NDREG=21) 

CHARACTER FORM*80,TITLE*64,NTC*2 

DIMENSION NRSUB(NDREG),ELFR(NDREG),ZET(NDREG) 

DIMENSION ALPI (NDREG) , BET I (NDREG) ,DELT (NDREG) 

DIMENSION NSG(NDREG),ZETG(NDREG) 

NT=LEN(TITLE) 

DO 80 1=1, NT 
IT-NT+1-I 

IF(TITLE(IT;IT).GT.' ')GO TO 81 

80 CONTINUE 

81 WRITE(NTC,'(I2)')IT 
URITEdFILE *)'&INPUTS' 

FORM='(4X,A8,2X,A1 ,A'//NTC//' ,A1 )' 

URITEdFILE, FORM)'TITLE TITLE,"" 


URITEd FILE, DMFACE *',IFACE,' ISIUN 


URITEd FILE, DMGROT =',IGROT,'NOI 


URITEdFILE, 2)'R0 =',R0,'EL 

URITEdFILE, 2)'RPM =',RPM,'RPM0 

URITEd FILE, 2)'PLEG =' ,PLEG, 'PRIG 

URITEd FILE, 2)'VISC =' ,VISC, 'DENS 

URITEdFILE, 2)'EMA *' ,EMA, 'ENA 

URITEdFILE, 2)'EMB =',EMB,'ENB 

URITEdFILE, 2)' HTAP *', HTAP, 'HBRL 

URITEdFILE, 2)'TOLH »',TOLH,'TOLV 

URITEdFILE, 3)' IHOME =',IHOME,'NITH 
NREG5=MIN(NREG,5) 

URITE(IFILE,4)'NREG =' ,NREG, ' NRSUB 


■', ISIUN 
»',NOI,'IFLOU 


*',EL,'C =' 

=',RPM0,'RPMD 
=',PRIG,'FZD 
=',DENS 
=',ENA 
*',ENB 
=', HBRL 
=',TOLV,'DUT 
I *',NITH,'NITV 


,C 

=' ,RPMD 
=',FZD 


URITE(IFILE,4)'NREG =' ,NREG, 'NRSUB =',(NRSUB( 
IF(NREG.GT.5)URITE(IFILE,6)(NRSUB(I),I=6,NREG) 
URITEd FI LE,5)'ELFR =',(ELFR(I),I=1,NREG5) 
IF(NREG.GT.5)WRITE(IFILE,7)(ELFR(I ), I=6,NREG) 
URITEd FI LE,5)'ZET =',(ZET(I),I=1,NREG5) 

IF(NREG.GT.5)URITE(IFILE,7)(ZET(I ), I=6,NREG) 
URITEd FILE, 5) 'ALPI =', (ALPI (I), 1=1, NREG5) 
IF(NREG.GT.5)URITE(IFILE,7)(ALPI(I),I=6,NREG) 
URITEd FI LE,5)'BETI *',(BETI(I),I=1,NREG5) 
IF(NREG.GT.5)URITE(IFILE,7)(BETI(I ), I=6,NREG) 
WRITE(IFILE,5)'DELT =',(DELT(I),I=1,NREG5) 
IF(NREG.GT.5)URITE(IFILE,7)(DELT(I).I=6,NREG) 
URITEd FILE,8)'NSG =',(NSG(I),I=1,NREG5) 

IF(NREG.GT.5)URITE(IFILE,9)(NSG(I),I=6,NREG) 
WRITE(IFILE,5)'ZETG =',(ZETG(I),I=1,NREG5) 
IF(NREG.GT.5)URITE(IFILE,7)(ZETG(I),I=6,NREG) 

1 FORMAT (4X,A8, I3,T30,A8, 13 ,T55, A8, 13) 

2 FORMAT(4X A8 1P,E12.4.T30,A8,E12.4.T55,A8,E12.4) 

3 FORMAT(4X,A8,I3,T30,A8,I4,T55,A8,I4) 

4 FORMAT(4X,A8,I3,T30,A8,5I4) 

5 FORMAT(4X,A8,1P,5E12.4) 

6 FORMAT(37X,5I4) 

7 FORMAT(12X,1P,5E12.4) 

8 FORMAT(4X,A8,5I4) 

9 FORMAT(12X,5I4) 

URITEdFILE,*)'/' 

RETURN 

END 


, ( NRSUB( I),I=1,NREG5) 



NASA/CR— 2003-212360 


04/14/1995 13:06 


Filename: SPIRAL I .FOR 


SUBROUT 1 HE TSEAL ( TOLV.N I TV, NOI , I FACE , ID I R , I GROT1 , NREG , NRSUB , 

+RE,REC,P1R,PIN,OMT,OMOT,OMDT,DUT,ZET,ALP, BET, DELT, ENGP, ZETG, 
+RGjzTG,H,U,V,P,TAU,CK,CB,AM,FLO,TOR,W, I AMASS, ITER, IER) 

C TURBULENT SEAL COMPUTATION SUBROUTINE 
C CALLED BY MAIN 
C CALLS VI SOLV, TORO, FORCE , KBCAL 
C FLAG DEFINITIONS: 

c NOI 9 1 NEGLECT AXIAL CONVECTIVE INERTIAL TERMS 

C I FACE = 1 FACE SEAL 

C 0 CYLINDRICAL SEAL 

C IDIR 9 1 AXIAL FLOW IS KNOWN POSITIVE 

C -1 AXIAL FLOW IS KNOWN NEGATIVE 

C I GROT = 1 GROOVES ROTATE 

C 0 GROOVES STATIONARY Mll> „ 

C -1 NO GROOVES (SET BY THIS SUB AND PASSED TO SUPPORTING SUBS) 

C I AMASS 9 1 CK, CB, AND AM ARE STIFFNESS, DAMPING AND MASS AT 0 FREQUENCY 

C 0 CK AND AM ARE DAMPING AND STIFFNESS AT DISTURBANCE FREQUENCY 

C I AMASS AND IER (ERROR COOE ) ARE OUTPUT 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

PARAMETER (NDZ=201 ,NDREG*21 ) 

DIMENSION NRSUB(NDREG),RG(NDZ,NDREG),ZET(NDREG),ZTG(NDZ,NDREG) 

DIMENSION ENGP(NDREG),ZETG(NDREG> 

DIMENSION ALP( NDREG) , BET ( NDREG) , SBET ( NDREG ) , CBET ( NDREG) , 
+DELT(NDREG),IGROT(NDREG),UHG(NDZ,NDREG),VHG(NDZ,NDREG) 

DIMENSION H(NDZ,NDREG),U(NDZ.NDREG),V(NDZ, NDREG), P(NDZ, NDREG) 

DIMENSION TAU(NDZ, NDREG), CK(4,4),CB(4,4),AM(4,4), 

+TMP(4,4) 

IER=0 

PI=4.D0*ATAN(1 .DO) 

C ADJUST FLAG AND GET SIN AND COS FOR SPIRAL GROOVE REGION 
DO 5 K=1 ,NREG . . „„ 

I F(ALP(K).LT.1 .D-8.0R.ABS(1 .DO-ALP(K) ) . LT. 1 .D-8.0R. 

+ ABS(BET(K)).LT.1.D-8)THEN 
C IGR0T=-1 SIGNIFIES NO GROOVES 
IGR0T(K)*-1 
ELSE 

IGR0T(K)=IGR0T1 

SBET(K)=SIN(BET(K)*PI/180.D0) 

CBET(K)=COS(BET(K)*PI/180.D0) 

END IF 
5 CONTINUE 

C HOME IN ON INITIAL VELOCITY AND GET VELOCITY AND PRESSURE DISTRIBUTIONS 
VI=ABS(PIN) 

VI=MAX(VI ,1 .D-3) 

CALL OUTSCR(' FIRST ORDER SOLUTION', 0) 

CALL VI SOLV( TOLV.N I TV, VI , NOI , I FACE , IDIR, NREG, NRSUB, 
+RE,REC,P1R,PIN,OMT,OMOT,DUT ,ZET, 

+ I GROT , ALP , SBET , CBET , DELT , ENGP, ZETG , UHG , VHG, 

+RG,ZTG,H,U,V,P, ITER, IER) 

I F( IER.NE . 0)RETURN 

C CALCULATE DIMENSIONLESS SHEAR STRESS AND FLOW AND TORQUE PARAMETERS 
FLO=RG( 1 , 1 )*H( 1 , 1 )*V( 1 , 1 )*2.D0*PI 

CALL TORQ(NOI, IFACE, IDIR, RE, REC.PIR.OMT, NREG, NRSUB, 

+IGROT.ALP, SBET, CBET, DELT, ENGP.UHG, VHG, 

+RG,ZTG,H,U,V,TAU,TOR> 

C CALCULATE LOAD UNDER FACE SEAL 

I F ( I FACE .EQ.DCALL FORCE (NREG , NRSUB , RG , ZTG , P , W) 

C CALCULATE 0 FREQUENCY STIFFNESS, CK, AND DAMPING 
CALL OUTSCR(' SECOND ORDER SOLUTION', 0) 

CALL KBCAL(NOI , I FACE, IDIR, NREG, NRSUB, 

+RE,REC,P1R,OMT,O.DO,DUT,ZET, 

+ I GROT , ALP , SBET , CBET , DELT , ENGP , ZETG , UHG, VHG , 

+RG,ZTG,H,U,V,CK,CB,IER) 

IFOER.NE.O )RETURN 

C IF NON ( 0 B VALUE T QF G D i STURBANc| M FREQ . AM AND CB WILL CONTAIN STIFFNESS AND 


C DAMPING AT DISTURBANCE FREQ. 

IAMASS=0 

CALL KBCAL(NOI , I FACE, IDIR, NREG, NRSUB, 

+ RE,REC,P1R,0MT,0MDT,DUT ,ZET, 

+ I GROT, ALP, SBET, CBET, DELT, ENGP, ZETG, UHG, VHG, 

+ RG,ZTG,H,U,V,AM,CB, IER) 

I F( IER.NE. 0)RETURN 
ELSE 

C OTHERWISE AM WILL CONTAIN MASS MATRIX AND CB WILL BE 0 FREQ DAMPING 
I AMASS 9 1 

CALL KBCAL(NOI , I FACE, ID I R, NREG, NRSUB, 

+ RE, REC.PIR.OMT , 1 .DO.DUT ,ZET , 

+ IGROT,ALP,SBET,CBET,DELT, ENGP, ZETG, UHG, VHG, 

+ RG,ZTG,H,U,V,AM,TMP,IER) 

I F ( I ER . NE . 0)RETURN 
DO 8 I=1,4-IFACE 
DO 8 J=1,4-IFACE 
8 AM(1,J) 9 CK(I,J)-AM(I,J) 

END IF 


END 
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FUNCTION FLMSHP(X) 

C THIS IS THE USER DEFINED FILM SHAPE FUNCTION 

C X IS THE DISTANCE FROM THE CENTER OF THE SEAL DIVIDED BY THE SEALING 
C LENGTH , L . -.5 <= X <= .5 

C FOR A SHAFT SEAL X = S/(2*L/D) (S IS ZTG IN COOE) 

C FOR A FACE SEAL X=(S+L/D-1)/(2*L/D) 

C FLMSHP IS THE SHAPE OF THE FILM (DIMENSIONAL) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/BFSHP/HTAP.HBRL 

FLMSHP=-HTAP*X+HBRL*(1.D0-(2.D0*X)**2) 

RETURN 

END 


FUNCTION FA(RE,H) 

C USER DEFINABLE FRICTION FACTOR FOR MOVING SURFACE 
C CALLED BY RECAL,TORQ,PHIPSI 
C COMMON BLOCK PASSED FROM MAIN PROGRAM 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
COMMON/BFAFB/EMA,ENA,EMB,ENB 
FA=MAX( 24 .DO/RE , ENA*RE**EMA) 

C H IS NOT USED NOW BUT MAY BE IN FUTURE FOR TREATING ROUGHNESS 
H1=H 
RETURN 
END 
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FUNCTION FB(RE,H) 

C USER DEFINABLE FRICTION FACTOR FOR STATIONARY SURFACE 
C CALLED BY RECAL,TORQ,PHIPSI 
C COMMON BLOCK PASSED FROM MAIN PROGRAM 
IMPLICIT DOUBLE PRECISION <A-H,0-Z) 
COMMON/BFAFB/EMA,ENA,EMB,ENB 
FB=MAX( 24 .DO/RE , ENB*RE**EMB ) 

C H IS NOT USED NOW BUT MAY BE IN FUTURE FOR TREATING ROUGHNESS 
H1=H 
RETURN 
END 


FUNCTION DELTP(RE,V,H,HSTEP,ZET) 

C USER DEFINABLE FUNCTION FOR COMPUTING LOSS COEFFICIENTS 

C PRESSURE CHANGE (DOWNSREAM -UPSTREAM) DUE TO SUDDEN CHANGE IN CROSS SECTION 
C HSTEP * STEP HEIGHT (H UPSTREAM-H DOWNSTREAM) 

C CALLED BY UVPCAL 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

IF(ABS(HSTEP).LT.1 .D-8)THEN 
DELTP=O.DO 
GO TO 99 

ELSE IF(HSTEP.LT.O.DO)THEN 
C COMPUTE LOSS COEFFICIENT FOR EXPANSION 
ZET1 =< 1 .DO- H/CH+HSTEP) )**2 
ELSE 

C USE INPUT LOSS COEFFICIENT FOR CONTRACTION (ZET) OR COMPUTE IT FROM RE 
RE1=RE*H*ABS(V) 

ZET1=ZET 
END IF 

DELTP*- ( 1 .D0+ZET1 ) 

I F( HSTEP . LT . 1 .D8)DELTP=DELTP+( H/ ( H+HSTEP ) )**2 
DELTP“DELTP*V*V 
99 RETURN 
END 


00 
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FUNCTION DIRFCN(RE,REC,P1R,0MT,VC0N, IFACE, 1D,R,H,U,DV) 
C CALCULATES PRIMARY FLOW DERIVATIVES FOR TAN VEL. U (ID=1) 

C OR PRESSURE P (ID=2) 

C CALLED BY UVPCAL 

n PA I I Q DUTDQT 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

V=VC0N/R/H 

CALL PHIPSI(IFACE,RE, REC,0MT,R,H,U,V f PHI ,PSI) 

I F ( I D . EQ . 1 )D I RFCN*-PH I /REC/V 
IF(ID.EQ.2)DIRFCN=-P1R*(PSI+REC*V*DV) 

RETURN 

END 
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SUBROUTINE RZGRIDCIFACE.ELT.NREG.NRSUB.ELFR.RG.ZTG) 

C GENERATES R (RZ) AND Z (ZTG) GRIDS 

C TO LOCATE STARTING POINT FROM Z ORIGIN ADD INSIDE RADIUS 
C FOR FACE SEAL (IFACE=1), SUBTRACT L/D FOR SHAFT SEAL 
C R = Z FOR FACE SEAL AND R=1 FOR SHAFT SEAL 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

PARAMETER (NDZ-201 , NDREG=21 ) 

DIMENSION NRSUB(NDREG),RG(NDZ,NDREG), ZTG(NDZ,NDREG),ELFR(NDREG) 
IF(IFACE.NE.1)IFACE=0 
ZTG(1,1)=-ELT 
RG(1,1)=1.D0 
I Ft I FACE. EQ. WHEN 
ZTG( 1,1 )=1 .DO-2 .D0*ELT 
RG(1,1) z ZTG(1,1) 

END IF 

DO 102 KK S 1,NREG 
IF(KK.GT.1)THEN 
ZTG( 1 , KK)=ZTG(NRS+1 ,KK- 1 ) 

RG(1,KK)*RG(NRS+1,KK-1) 

END IF 

NRS=NRSUB(KK) 

DZF=2.D0*ELFR(KK)*ELT/NRS 
DO 102 JJ*1,NRS 

ZTG(JJ+1,KK)=ZTG(JJ,KK)+DZF 
RG( JJ+1 ,KK)»1 .DO 

I F( I FACE . EQ . 1 >RG< J J+1 , KK)=ZTG< J J+1 , KK) 

102 CONTINUE 
RETURN 
END 
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SUBROUTINE VISOLV(TOLV,NITV, VI, NOI, I FACE, IDIR, NREG, NRSUB, 
+RE,REC,P1R,PIN,OMT,OMOT,DUT,ZET, 

+ I GROT , ALP , SBET , CBET , BELT , ENGP , ZETG , UHG , VHG , 

+RG,ZTG,H,U,V,P, ITER, IER) 

C SOLVES FOR INLET VELOCITY VI USING NEWTONS METHOD 
C ON INPUT VI IS INITIAL GUESS. IER=2 IF NOT CONVERGED. 

C CALLED BY TSEAL 
r TAI I ^ UVPCAL 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

PARAMETER <NDZ=201,NDREG=21) 

DIMENSION NRSUB( NDREG), RG(NDZ, NDREG), ZET(NDREG),ZTG(NDZ, NDREG) 
DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ,NDREG),P(NDZ,NDREG) 
DIMENSION ALP( NDREG) , SBET(NDREG) , CBET (NDREG) ,DELT(NDREG ) , 

+ 1 GROT ( NDREG ) , UHG( NDZ, NDREG) , VHG< NDZ , NDREG > 

DIMENSION ENGP(NDREG),ZETG(NDREG) 

DO 5 J=1,NITV 
ITER=J 

CALL UVPCAL ( VI +DUT , NOI , I FACE , ID I R , NREG , NRSUB , 

+ RE,REC,P1R,PIN,OMT,OMOT,DUT,ZET , 

+ I GROT , 6 , ALP , SBET .CBET , DELT , ENGP , ZETG, UHG , VHG , 

+ RG,ZTG,H,U,V,P,F2, IER) 

CALL UVPCAL (VI , NOI , I FACE , I D I R , NREG , NRSUB , 

+ RE,REC,P1R,PIN,OMT,OMOT,DUT,ZET, 

+ I GROT , 1 , ALP , SBET . CBET , DELT , ENGP , ZETG , UHG , VHG , 

+ RG,ZTG,H,U,V.P,F1,IER) 

DV=F1*DUT/(F2-F1) 

VI=VI-DV 
VI 1=ABS(VI )+TOLV 
IF(ABS(DV)/VI1.LT.T0LV)G0 TO 6 

5 CONTINUE 
IF(NITV.GT.1)IER S 2 

6 CALL UVPCAL(VI, NOI, IFACE.IDIR, NREG, NRSUB, 
+RE,REC,P1R,PIN,OMT,OMOT,DUT,ZET, 

+IGROT.1, ALP, SBET. CBET, DELT, ENGP, ZETG, UHG, VHG, 
+RG,ZTG,H,U,V,P,F1, IER) 

RETURN 

END 
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SUBROUTINE UVPCAL(VI ,NOI , IFACE, IDIR, NREG, NRSUB, 
+RE,REC,P1R,PIN,OMT ,OMOT ,DUT,ZET , 

+IGROT,IUHG,ALP, SBET, CBET, DELT, ENGP, ZETG,UHG,VHG, 
+RG,ZTG,H,U,V,P,PEXIT,IER) 

C GENERATES PRESSURE (P) AND VELOCITIES (U,V) BASED ON INITIAL VALUE OF 

C V (VI) 

C CALLED BY VISOLV 

C CALLS UVPIN WHEN REC>0 OR UVPNOI WHEN REC = 0 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

PARAMETER (NDZ=201 ,NDREG=21 ) „ 

DIMENSION NRSUB(NDREG),RG(NDZ, NDREG), ZET(NDREG),ZTG(NDZ, NDREG) 
DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ, NDREG), P(NDZ,NDREG) 
DIMENSION ALP(NDREG),SBET(NDREG),CBET(NDREG),DELT(NDREG), 
+IGROT(NDREG),UHG(NDZ,NDREG),VHG(NDZ, NDREG) 

DIMENSION ENGP(NDREG),ZETG(NDREG) 

I F( NOI .EQ.DTHEN 

CALL UVPNOI (VI , I FACE, IDIR, NREG.NRSUB, 

+ RE,REC,P1R,PIN,OMT,DUT, 

+ I GROT, IUHG, ALP, SBET, CBET, DELT, ENGP, ZETG, UHG, VHG, 

+ RG,ZTG,H,U,V,P,PEXIT,IER) 

ELSE 

CALL UVPI N ( VI , I FACE , I D I R , NREG , NRSUB , 

+ RE,REC,P1R,PIN,OMT,OMOT,DUT,ZET , 

+ I GROT, IUHG,ALP,SBET, CBET, DELT, ENGP.ZETG.UHG, VHG, 

+ RG,ZTG,H,U,V,P f PEXIT.IER) 

END IF 

RETURN 

END 
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SUBROUT I NE UVPNOI ( VI , I FACE , I D I R , NREG , NRSUB, 

+RE,REC,P1R,PIN,0MT,DUT, 

+IGR0T, IUHG, ALP, SBET,CBET,DELT,ENGP,ZETG,UHG,VHG, 

+RG,ZTG,H,U,V,P,PEXIT, IER) 

C GENERATES PRESSURE (P) AND VELOCITIES (U,V) BASED ON INITIAL VALUE OF 
C V (VI) WITHOUT INERTIA EFFECTS 
C CALLED UVPCAL 

C CALLS DIRFCN.USOLV (NO GROOVES) OR PHIPSG (FOR SPIRAL GROOVES) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

PARAMETER (NDZ=201,NDREG=21) 

C COMMON BLOCK USED LOCALLY IN THIS ROUTINE 

COMMON/BGRLCL/UHL(NDZ , NDREG ),VHL(NDZ, NDREG ) 

DIMENSION NRSUB(NDREG), RG(NDZ, NDREG), ZTG(NDZ, NDREG) 

DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ,NDREG),P(NDZ, NDREG) 

DIMENSION ALP(NDREG) ,SBET (NDREG) ,CBET (NDREG) ,DELT(NDREG), 

+1 GROT (NDREG) ,UHG(NDZ, NDREG), VHG(NDZ, NDREG) 

DIMENSION ENGP(NDREG),ZETG(NDREG) 

C SET UP STARTING CONDITIONS AND LOOPING PARAMETERS FOR PRIMARY FLOW SOLUTION 
KST=1 
KEN=NREG 

I F( IDIR.EQ. - 1 )TKEN 
KST=NREG 
KEN=1 
END IF 

DO 30 K-KST,KEN,IDIR 
JST=1 

JEN=NRSUB(K) 

I F( IDIR.EQ. -1)THEN 
JST=JEN+1 
JEN=2 
END IF 

IF(K.EQ.KST)THEN 

VCON=VI *RG( JST , KST )*H ( JST , KST )* I D I R 

V(JST,KST)=VI*IDIR 

P(JST,KST)=PIN 

ELSE 

V(JST,K)=VCON/RG(JST,K)/H(JST,K) 

P(JST,K)=P( J1 ,K1 ) 

END IF 

IF(IGROT(K).EQ.-1)THEN 

CALL USOLV(RE,OMT,RG( JST,K),H( JST,K),U(JST,K), 

+ V(JST,K),DUT, IUHG, IER) 

ELSE 

CALL PHIPSG(1 , 1 FACE, IUHG, RE, REC,0MT,RG( JST, K),H( JST, K), 

+ U(JST,K),V(JST,K),UHL(JST,K),VHL(JST,K), 

+ ALP(K),SBET(K),CBET(K),DELT(K),ENGP(K),ZETG(K), 

+ IGROT(K),PHI ,PSI , IER) 

END IF 

IF(IER.NE.O)RETURN 

K1=< 

DO 30 J=JST,JEN,IDIR 
J1*J+IDIR 

RB=.5D0*(RG(J1,K)+RG(J,K)) 

HB=.5D0*(H( J1 ,K)+H( J,K)) 

DX=ZTG(J1,K)*ZTG( J,K) 

V(J1,K)=VCON/RG(J1,K)/H(J1,K) 

IF(IGROT(K).EQ.-1)THEN 

CALL US0LV(RE,0MT,RG(J1,K),H(J1,K),U(J1,K),V(J1,K), 

+ DUT, IUHG, IER) 

I F ( 1 ER . NE . 0 )RETURN 

P(J1,K)=P(J,K)+DX*DIRFCN(RE,REC,P1R,0MT,VC0N,IFACE,2, 

+ RB,HB,.5D0*(U(J1,K)+U(J,K)),0.D0) 

ELSE 

VB=VCON/RB/HB 

CALL PHIPSG(1 ,IFACE, IUHG,RE,REC,OMT,RG(J1 ,K),H(J1,K), 

+ U(J1,K),V(J1,K),UHL(J1,K),VHL(J1,K), 



+ ALP(K),SBET(K),CBET(K),DELT(K),ENGP(K),ZETG(K), 

+ IGROT(K),PHI,PSI,IER) 

IF(IER.NE.O)RETURN 

CALL PH I PSG(0, 1 FACE , IUHG, RE , REC.OMT , RB, HB , 

+ .5D0*(U(J1,K)+U(J,K)),VB,UHG(J,K),VHG(J,K>, 

+ ALP(K),SBET(K),CBET(K),DELT(K),ENGP(K),ZETG(K), 

+ IGROT(K),PHI,PSI,IER) 

P(J1,K)=P(J,K)-DX*P1R*PSI 
IF(IER.NE.O)RETURN 
END IF 
30 CONTINUE 

PEXIT=P(J1,K1) 

RETURN 


END 
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SUBROUTINE USOLV(RE,OMT,R,H,U,V,DUT, HAST, IER) 

IMPLICIT DOUBLE PRECISION <A-H f 0-Z) 

C SOLVES FOR EQUILIBRIUM TANGENTIAL VELOCITY WHEN THERE ARE NO GROOVES 
C CALLED BY UVPNOI 
C CALLS PHIPSI 

IF(ABS(OMT).LT.DUT)THEN 
U=0.D0 
RETURN 
END IF 

DU=DUT*0MT 

TOL=100.D0*ABS(DU) 

I F ( I LAST . EQ . 0 )U= . 5D0*0MT 
DO 5 1=1,30 

CALL PHIPSI(0,RE,O.DO,OMT,R,H,U,V,PHI ,PSI ) 

CALL PHIPSI(0,RE,O.DO,OMT,R,H,U+DU,V,DPHI ,PSI ) 
DLT=PHI*DU/(DPHI-PHI) 

U=U-DLT 

IF(ABS(DLT).LT.TOL)GO TO 6 

5 CONTINUE 
IER=2 

6 RETURN 
END 


SUBROUT I NE UVPI N( VI , I FACE , I D I R, NREG , NRSUB , 

+RE,REC,P1R,PIN, OMT, OMOT ,DUT,ZET , 

+IGROT, IUHG,ALP,SBET,CBET,DELT,ENGP,ZETG,UHG,VHG, 
+RG,ZTG,H,U,V,P,PEXIT,IER) 

C GENERATES PRESSURE (P) AND VELOCITIES (U,V) BASED ON INITIAL VALUE OF 
C V (VI) WHEN INERTIA IS PRESENT 
C CALLED BY UVPCAL 

C CALLS DIRFCN,DELTP,PHIPSG (FOR SPIRAL GROOVES) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

PARAMETER (NDZ=201,NDREG=21) 

C COMMON BLOCK USED LOCALLY IN THIS ROUTINE 

COMMON/BGRLCL/UHL(NDZ,NDREG),VHL(NDZ,NDREG) 

DIMENSION NRSUB(NDREG),RG(NDZ,NDREG),ZET(NDREG),ZTG(NDZ,NDREG) 
DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ,NDREG),P(NDZ,NDREG) 

DIMENSION ALP(NDREG) , SBET(NDREG) , CBET(NDREG) ,DELT(NDREG) , 
+IGROT(NDREG),UHG(NDZ,NDREG),VHG(NDZ, NDREG) 

DIMENSION ENGP(NDREG),ZETG(NDREG) Bnl ,„, nu 

C SET UP STARTING CONDITIONS AND LOOPING PARAMETERS FOR PRIMARY FLOW SOLUTION 
KST=1 
KEN=NREG 

RV20=.5D0*REC*P1R 
IF(IDIR.EQ.-1)THEN 
KST=NREG 
KEN=1 
END IF 

DO 30 K=KST,KEN,IDIR 
JST*1 

JEN=NRSUB(K) 

IF(IDIR.EQ.-1)THEN 
JST=JEN+1 
JEN=2 
END IF 

IF(K.EQ.KST)THEN 

VCON=VI*RG(JST,KST)*H(JST,KST)*IDIR 

V(JST,KST)=VI*IDIR 

U( JST , KST )=OMOT*RG( JST , KST) 

P(JST,KST)=PIN+ 

+ RV20*DELTP( RE , V( JST , KST ) , H( JST , KST ),1.D10,ZET( KST ) ) 

ELSE 

V( JST,K)=VCON/RG( JST ,K)/H( JST,K) 

U(JST,K)=U(J1,K1) 

P(JST,K)=P(J1,K1)+ 

+ RV20*DELTP(RE,V(JST,K),H( JST,K),H( J1,K1)-H( JST ,K),ZET(K)) 

END IF 
K1=K 

DO 30 J=JST,JEN, IDIR 
J 1 = J+ I D I R 

RB=.5D0*(RG(J1 ,K)+RG(J,K)) 

HB=.5D0*(H(J1,K)+H(J,K)) 

DX=ZTG(J1,K)-ZTG(J,K) 

V( J1 ,K)=VCON/RG( J1 ,K)/H( J1 ,K) 

IF(IGROT(K).EQ.-1)THEN 

U1=DIRFCN(RE,REC,P1R,0MT,VC0N,IFACE,1, 

+ RB,HB,U(J,K),0.D0) 

DU1=(DIRFCN(RE,REC,P1R,0MT,VC0N,IFACE,1 , 

+ RB,HB,U(J,K)+DUT,0.D0)-U1)/DUT 

U( J1 ,K)=U( J ,K)+DX*Ul/( 1 .DO- .5D0*DX*DU1 ) 

P( J1 ,K)=P(J f K)+OX*DIRFCN(RE,REC,P1R,OMT,VCON, IFACE.2, 

+ RB,fiB,.5D0*(U(J1,K)+U(J,K)),(V(J1,K)-V(J,K))/DX) 

ELSE 

VB=VCON/RB/HB 

CALL PHIPSG(0,IFACE,IUHG,RE f REC,OMT,RB,HB, 

+ U( J, K),VB,UHL(J,K),VHL(J,K), 

+ ALP(K), SBET(K),CBET(K),DELT(K),ENGP(K),ZETG(K), 

+ IGROT(K),PHI,PSI,IER) 

IF(IER.NE.O)RETURN 
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U1=-PHI/REC/V(J,K) 

CALL PHIPSG<0,IFACE f -1,RE,REC,OMT,RB f HB, 

+ U( J,K)+DUT,VB,UHL(J,K),VHL(J,K), 

+ ALP(K),SBET(K),CBET(IO < DELT(K),ENGP(K),ZETG<K), 

+ IGROT(K),PHI,PSI,IER) 

IF(IER.NE.O)RETURN 

DU1 =< -PH I /REC/VB-U1 )/0UT 

U(J1,K)=U(J,K)+DX*U1/(1.D0-.5D0*DX*DU1) 

CALL PHIPSG(0, IFACE,IUHG,RE,REC,OMT,RB,HB, 

+ .5D0*<U(J1,K)+U<J,K)) f VB,UHG(J,K),VHG(J,K) f 

+ ALP(K),SBET<K),CBET(K),DELT<K),ENGP(K),ZETG<K), 

+ IGROT(K),PHI ,PSI , IER) 

P<J1,K)=P(J t K>-DX*P1R*<PSI+REC*VB*(V<J1,K)-V(J,K))/DX) 
IF(IER.NE.O)RETURN 
END IF 

30 CONTINUE 

PEXIT=P(J1 f K1) 

RETURN 

END 



SUBROUTINE RECAL(REL,RE,FBAR,NITV,TOLV, IER) 

C USES NEWTON ITERATION TO GET CHARACTERISTIC REYNOLDS NUMBER FOR TURBULENT 
C POISEUILLE FLOW WITH UNIFORM CLEARANCE 
C REL = LAMINAR REYNOLDS NUMBER (INPUT) 

C RE = REYNOLDS NUMBER 
C IER - ERROR CODE, 0 IF OK 
C CALLS FA,FB 
C CALLED BY MAIN PROGRAM 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

IER=0 

RE=REL 

DRE=1 .D-6*RE 
DO 5 1=1 ,NITV 
11=1 

FBAR=(FA(RE,1.D0)+FB(RE,1.D0))/2.D0 

DFBAR=((FA(RE+DRE.1.D0)+FB(RE+DRE,1.D0))/2.D0-FBAR)/DRE 

DELT=(RE*RE*FBAR-24.D0*REL)/(RE*RE*DFBAR+2.D0*RE*FBAR) 

RE=RE-DELT 

IF(ABS(DELT/RE).LT.TOLV)GO TO 6 
5 CONTINUE 


END 
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SUBROUT I NE FORCE ( NREG , NRSUB , RG , ZTG, P, W) 

C COMPUTES DIMENSIONLESS LOAD, W, FROM PRIMARY PRESSURE DISTRIBUTION 
C ONLY MEANINGFUL FOR FACE SEAL 
C CALLED BY TSEAL 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

PARAMETER (NDZ=201,NDREG=21) 

DIMENSION NRSUB(NDREG),RG(NDZ,NDREG), ZTG(NDZ, NDREG), P(NDZ, NDREG) 
W=0.D0 

DO 5 K-1,NREG 
NRS-NRSUB(K) 

DO 5 J=1 ,NRS 

J1=J+1 , „ 

5 U=W+(P(J,K)+P( J1 ,K))*(RG(J,K)+RG( J1 ,K))*(ZTG(J1, K)-ZTG( J,K))/4.D0 

W=W*8.D0*ATAN(1.D0) 

RETURN 

END 
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SUBROUTINE TORQ(NOI / IFACE,IDIR,RE,REC f P1R,OMT, NREG, NRSUB, 

+1 GROT, ALP, SBET,CBET,DELT,ENGP,UHG,VHG, 

+RG,ZTG,H,U,V,TAU,TOR) 

C CALCULATES SHEAR STRESS ON MOVING SURFACE AND TORQUE INTEGRAL 
C SHEAR STRESSES AT ARE AT HALF GRID POINTS 
C CALLED BY TSEAL 
C CALLS FA AND FB 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

PARAMETER (NDZ=201,NDREG=21> 

DIMENSION NRSUB(NDREG),RG(NDZ,NDREG),ZTG(NDZ,NDREG),TAU(NDZ,NDREG) 
DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ,NDREG) 

DIMENSION ALP(NDREG) ,DELT(NDREG) , SBET (NDREG) .CBET(NDREG) , 
+IGROT(NDREG),UHG(NDZ,NDREG),VHG(NDZ, NDREG) 

DIMENSION ENGP( NDREG) 

C REDUCE COUETTE PART OF SHEAR STRESS BY 3 FOR LAMINAR FLOW 
C OR BY 1.2 AS SUGGESTED BY HIRS FOR TURBULENT FLOW 
DATA DLAM.DHIRS/3.D0.1 .200/ 

VCON=V< 1 , 1 )*RG( 1 , 1 )*H< 1 , 1 ) 

TOR-0. DO 

C NEED TO SET INTEGRATION DIRECTION TO PROPERLY GET UHG AND VHG AT HALF GRID 
KST=1 
KEN=NREG 

IF(IDIR.EQ.-1)THEN 
KST=NREG 
KEN=1 
END IF 

DO 4 K-KST,KEN,IDIR 
JST=1 

JEN=NRSUB(K) 

IF(IDIR.EQ.-1)THEN 

JST=JEN+1 

JEN=2 

ENDIF 

DO 4 J=JST, JEN.IDIR 
J1=J+IDIR 

HBAR=.5D0*(H( J,K)+H(J1,K)) 

UBAR=.5D0*(U(J,K)+U(J1,K)) 

RBAR=.5D0*(RG(J,K)+RG(J1,K)) 

VBAR-VCON/RBAR/HBAR 

DPCOR=O.DO 

IF(IGROT(K).EQ.-1)THEN 
RA=RE*HBAR*SQRT ( (UBAR - RBAR*OMT )**2+VBAR**2 ) 

RAFA=RA*FA(RA,HBAR) 

C REDUCE SHEAR STRESS BY FACTOR OF 3 FOR LAMINAR FLOW 
DLM-DHIRS 

IF(ABS(RAFA-24.D0).LT.1.D-10)DLM=DLAM 

TAU(J,K)=P1R*RAFA*(UBAR-RBAR*OMT)/HBAR/DLM 

ELSE 

HR=HBAR-ALP(K)*DELT(K) 

HG=HR+DELT(K) 

UHR=(UBAR*HBAR-UHG(J,K)*ALP(K))/(1 .DO-ALP(K)) 
VHR=(VBAR*HBAR-VHG(J,K)*ALP<K))/(1 .DO-ALP(K)) 

RAG=RE*SQRT< (UHG( J , K) -RBAR*0MT*HG)**2+VHG( J ,0**2) 

RAR=RE*SQRT ( (UHR-RBAR*OMT*HR >**2+VHR**2 ) 

RAFAG=RAG*FA(RAG, HG) 

RAFAR=RAR*FA(RAR,HR) 

RBG=RE*SQRT(UHG(J,K)**2+VHG(J,K)**2) 

R8R=RE*SQRT (UHR**2+VHR**2 ) 

RBFBG=RBG*FB(RBG,HG) 

RBFBR=RBR*FB(RBR,HR) 

TAUGA=RAFAG*(UHG(J,K)-RBAR*0MT*HG)/HG**2 
TAUGB*-RBFBG*UHG( J , K)/HG**2 
TAURA«RAFAR*(UHR-RBAR*0MT*HR)/HR**2 
TAURB=-RBFBR*UHR/HR**2 
C SET COUETTE REDUCTION FACTOR 
DLMsDHIRS 
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IF(ABS<RAFAR-24.D0).LT.1.D-10.OR.ABS(RAFAG-24.D0).LT.1.D-10) 

+ DLH=DLAM 

C ATTEMPT TO SPLIT OFF COUETTE AND POISEUILLE PORTIONS OF SHEAR STRESS 
T AURC= ( T AURA+T AURB )/2 . DO/D LM 
T AUGC= < T AUGA+T AUGB )/2 . DO/DLM 
TAURP=( TAURA- TAURB )/2 . DO 
TAUGP=< TAUGA- TAUGB )/2 . DO 

C IF GROOVES ROTATE CORRECT FOR FORCES AT GROOVE EDGES 
IF(IGROT(K).EQ.1)THEN 
TAUGP=TAUGP*( 1 . DO-2 . D0*DELT <K)/HG ) 

C GET EFFECTS OF LOCAL INERTIA DROP, DPCOR 

QN=HBAR*(UBAR-RBAR*OMT)*SBET(K) 

IF(NOI.NE.1 )QN=QN-HBAR*VBAR*CBET (K) 

DPCOR=REC*ENGPOO/RBAR*DELTP(RE,QN/HG,HG,-DELT(K),O.DO)* 

+ SIGN(1.D0,SBET(K)*QN) 

C 

IF(NOI.NE.1)THEN 

TICOR=(U(J1,K)-U( J,K))/(ZTG(J1,K)-ZTG( J,K)) 

I F ( I FACE . EQ . 1 )T I COR=T I COR+UBAR/RBAR 
TAUGP=TAUGP-DELT<K)*REC*VBAR*TICOR 
END IF 
END IF 

TAU(J,K)=P1R*(ALP(K)*(TAUGC+TAUGP)+ 

+ (1 .DO-ALP(K))*(TAURP+TAURC)+DPCOR) 

END IF 

TOR=TOR+TAU( J,K)*RBAR**2*(ZTG( Jl ,K)-ZTG(J,K))*IDIR 
4 CONTINUE 

C MULTIPLY BY 2 PI AND CHANGE SIGN SO THAT TORQUE IS + WHEN IT OPPOSES MOTION 
TOR*- T0R*8 . DO*AT AN ( 1 .DO) 

RETURN 
END 
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SUBROUTINE PHIPSG(NOI, IFACE, IUHG,RE,REC,OMT,R,H / U,V,UHG,VHG, 
+ALP,S,C,DELT,ENGP,ZETG, IGROT.PHI ,PSI , IER) 

C GENERATES GLOBAL TURBULENCE FUNCTIONS FOR SPIRAL GROOVES 
C CALLED BY UVPNOI,UVPIN,DSOLV 
C CALLS PHIPSQ.MATINV 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION Q(4), B(4),INDEX<4,3), A(4,4),D(4,4),E(4) 

DATA A/16*0.D0/ 

DQ=1.D-6 

HR=H-ALP*DELT 

HG=HR+DELT 

UH=U*H 

VH*V*H 

IF(IUHG.EQ.O)THEN 

I F(NOI .EQ. 1 )UH=.5DO*H*R*OMT 

Q(1)=UH 

Q(2)*VH 

ELSE 

Q(1)4JHG 
Q(2)=VHG 
END IF 

Q(3)=(UH-Q<1 )*ALP)/( 1 .DO-ALP) 

Q(4)s(VH-Q(2)*ALP)/( 1 .DO-ALP) 

DO 30 L=1,30 

CALL PHIPSQ(RE,OMT ,R,HG,Q(1),Q(2),B(1),B(2)) 

CALL PHIPSQ(RE,OMT,R,HG,Q(1)+DQ,Q(2),A(1,1),A(2,1)) 

CALL PHIPSQ(RE,0MT,R,HG,QO),Q(2)+DQ,AO,2),A<2,2)) 

CALL PHIPSQ(RE,OMT,R,HR f Q(3),Q(4),B<3),B(4)) 

CALL PHIPSQ(RE,OMT,R,HR,Q(3)+DQ,Q(4),A(3,3),A(4,3)) 

CALL PHIPSQ(RE,OMT,R, HR,Q(3),Q(4)+DQ,A(3,4),A(4,4)) 

DO 8 K1=1,2 
X=2*(K1-1) 

DO 8 1*1,2 
DO 7 J=1 ,2 

7 A(I+K,J+K)=(A(I+K,J+K)-B(I+K))/DQ 
B(I+K)*-B(I+K) 

DO 8 J*1,2 

8 B( I+K)=B( I+K)+A( I+K, J+K)*Q( J+K) 
E(1)*C*B(1)+S*B(2)-C*B(3)-S*B(4) 

DO 9 J=1 ,4 

9 D(1 , J)*C*A(1, J)+S*A(2, J)-C*A(3,J)-S*A(4, J) 

D(2,1)=S 

D<2,2)*-C 

D(2,3)*-S 

D(2,4)*C 

E ( 2 ) *R*OMT*DE LT*S* I GROT 
IF(NOI.EQ.1)THEN 
E(3)*ALP*B(1)+(1.D0-ALP)*B(3) 

DO 10 J*1,4 

10 D(3,J)=ALP*A(1,J)+(1.D0-ALP)*A(3,J) 

C IF(IFACE.EQ.1)THEN 

C COR=REC*VH/R/H**2 

C D(3, 1 )=D(3, 1 )-ALP*COR 

C D(3,3)=D(3,3)-(1.D0-ALP)*C0R 

C END IF 

ELSE 

D(3,1)*ALP 
D(3,2)*0.D0 
D(3,3)=1. DO-ALP 
D(3,4)*0.D0 
E(3)=UH 
END IF 

D(4,1)*0.D0 
D(4,2)=ALP 
D(4,3)*0.D0 
D(4,4)=1 .DO-ALP 
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E(4)=VH 

CALL MATINV(D,E,DETER,4,1 , ID, 4, INDEX) 

I F( ID . NE . 1 )G0 TO 99 
I CNV=0 
EMX=-1 .020 
00 11 1=1,4 

EMX=MAX(EMX,ABS(E( I)-Q(l))) 

11 Q(I)=E(I) 

IF(EHX.LT.1.D-4)G0 TO 31 

30 CONTINUE 
GO TO 99 

31 CALL PHIPSQ(RE,0MT,R,HG,Q(1),Q(2),PHIG,PSIG) 

CALL PHIPS0(RE,0MT,R,HR,Q(3),Q(4), PHIR,PSIR) 

I F(N0I .EQ. 1 )U=(ALP*Q(1 )+< 1 .D0-ALP)*Q(3))/H 

I F< IUHG.NE . - 1 )THEN 
UHG=Q(1) 

VHG=Q(2) 

END IF 

E RR=C*PH I G+S*PS I G - C*PH I R - S*PS I R 
IF(ERR.GT.1.D-4)IER=4 
PHI=ALP*PHIG+(1 .DO-ALP)*PHIR 
PSI=ALP*PSIG+(1 .DO-ALP)*PSIR 
C ADD EFFECTS OF LOCAL INERTIA DROP 
C NORMAL FLOW 

QN=H*(U-R*OMT*IGROT)*S 
I F(N0I .NE. 1 )ON=QN-H*V*C 
C CONTRACTION LOSS COEFF, ZETG 

PDRP=DELTP(RE.QN/HG.HG,HR-HG,0.D0)+DELTP(RE,QN/HR,HR,HG-HR,2ETG) 

PH I =PH I - S I GN ( 1 .DO , S*QN )*REC*ENGP/R*PDRP 

PS I =PS I+S I GN ( 1 .DO , C*QN )*REC*ENGP/R*PDRP*ABS(C/S ) 

I F ( I FACE . NE . 1 )RETURN 
I F(NOI .NE. 1 )PH I =PH I +REC*U*V/R 
PSI*PSI* REC*U*U/R 
RETURN 
99 IER=4 
RETURN 
END 


SUBROUTINE PHIPSQ<RE,OMT,R,H,QT,QS,PHI,PSI) 

C GENERATES TURBULENCE FUNCTIONS PHI.PSI BASED ON FLOW RATHER THAN VELOCITY 
C EXCLUDES CENTRIFUGAL AND CORIOLIS TERMS FOR FACE SEAL 


C CALLS PHIPSI 
C CALLED BY PHIPSG 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

CALL PHIPSI(0,RE, O.DO,OMT, R,H,QT/H,QS/H,PHI,PSI ) 

RETURN 

END 
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SUBROUTINE PHIPSI(IFACE,RE,REC,OMT,R,H,U,V,PHI,PSI) 
C GENERATES TURBULENCE FUNCTIONS PHI.PSI 
C INCLUDES CENTRIFUGAL AND CORIOLIS TERMS FOR FACE SEAL 
C CALLS FA FB 

C CALLED BY DIRFCN,DSOLV,PHIPSQ 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
RA=RE*H*SQRT((U-R*0MT)**2+V*V) 

RB=RE*H*SQRT (U*U+V*V) 

CK2=RA*FA(RA,H) 

CK1 =CK2+RB*FB( RB , H ) 

PHI=(CK1*U-CK2*R*0MT)/H**2 

PSI=CK1*V/H**2 

I Ft I FACE . NE . 1 )RETURN 

PH I =PH I +REC*U*V/R 

PS I =PS I - REC*U*U/R 

RETURN 

END 
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SUBROUTINE DSOLV(NOI , I FACE,RE,REC1 ,P1R,0MT,0MDT ,DUT , 

+ALP, SBET , CBET , DELT , ENGP , ZETG, I GROT , 
+R1,Z1,H1,U1,V1,UHG,VHG,R,Z,H,U,V,Y,F0I,EMI,IER) 

C UPDATES Y(I,L.N) TO NEXT Z POSITION 

C GENERATES MATRICES FOR DISURBANCE EQUATIONS AT ONE VALUE OF Z 
C DISTURBANCE EQUATIONS ARE IN FORM {DY/DZ> = [AKY> + CB> 

C WHERE 1=1, 2,3 CORRESPONDS TO P,V,U DISTURBANCES 

C L=1,2,3 CORRESPONDS TO COMPLIMENTARY SOLUTION, TILT.RADIAL DISP RESP. WHEN N<3 
C N=1,2,3 CORRESPONDS TO EXP(I*(THETA+OMDT*T)),EXP(I*(THETA-OMOT*T)), 

C EXP(I*OMDT*T) RESP. 

C N=3 IS FOR AXIAL DISTURBANCE APPLIED TO FACE SEAL FOR WHICH CASE 
C L=1,2 CORRESPONDS TO COMP. SOL. AND AXIAL DISP RESP. 

C CALLS PHIPSI,PHIPSG,ESET, CHAT I N 
C CALLED BY KBCAL 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

C0MPLEX*16 A(3,3),B(3,3), TI(3),0I(3),UR,DETER,Y(3,3,3), 
+F0I(3,3),EMI<3,3),TMP,DFI 
DIMENSION E(3),DE(3),LABEL(3,3) ^ „ 

DATA TI,OI/<O.DO,1.DO),(O.DO,1.DO),<O.DO.O.DO), 

+ (O.DO f 1.DO),(O.DO,-1.DO),(O.DO,1.DO)/ 

IF(NOI?EQ.1)REC=O.DO 

UB=.5D0*(U+U1) 

HB=.5D0*(H+H1) 

RB=.5D0*(R+R1) 

VB=V*R*H/RB/HB 

ZB=.5D0*(Z+Z1) 

DZ=Z-Z1 

DZ2=.5D0*DZ 

DU=(U-U1)/DZ 

DV=(V-V1)/DZ 

DRH=(RB*(H-H1 )+HB*(R-R1 ) )/DZ 
REP=P1R*REC 

C CALCULATE TURBULENT FUNCTIONS AND THEIR DERIVATIVES 
I F( I GROT .EQ.-1 )THEN 

CALL PHIPSI(IFACE,RE,REC,OMT,RB,HB,UB,VB,PHI ,PSI) 

CALL PHIPSI(IFACE,RE,REC,OMT,RB,HB+OUT,UB,VB,PHIH,PSIH) 

CALL PHIPSI(IFACE,RE,REC,OMT,RB,HB,UB+DUT,VB,PHIU,PSIU) 

CALL PHIPSI(IFACE,RE,REC,OMT,RB,HB,UB,VB+DUT,PHIV,PSIV) 

ELSE 

UHGB=UHG 

VHGB=VHG 

CALL PHIPSG(0,IFACE, 1 ,RE,REC,OMT,RB,HB,UB,VB,UHGB,VHGB, 

+ ALP, SBET, CBET, DELT, ENGP, ZETG, IGROT, PHI ,PSI , IER) 

IFCIER NE 0)RETURN 

CALL PHIPSG(0, IFACE.-1 ,RE,REC,OMT,RB,HB+DUT,UB,VB,UHGB,VHGB, 

+ ALP, SBET, CBET, DELT, ENGP, ZETG, IGROT ,PHIH,PSIH, IER) 

IFCIER NE 0)RETURN 

CALL PH I PSGCO, I FACE ,-1,RE,REC, OMT , RB , HB , UB+DUT , VB , UHGB , VHGB , 

+ ALP, SBET, CBET, DELT, ENGP, ZETG, IGROT, PHIU.PSIU, IER) 

1F(IER.NE.0)RETURN 

CALL PHIPSG(0, 1 FACE, - 1 ,RE,REC,OMT,RB,HB,UB,VB+DUT, UHGB, VHGB, 

+ ALP, SBET, CBET, DELT, ENGP.ZETG, IGROT, PHIV,PSIV, IER) 

IF(IER.NE.O)RETURN 
END IF 

PHIH=(PHIH-PHI)/DUT 
PSIH=(PSIH-PSI)/DUT 
PHIU=(PHIU-PHI)/DUT 
PSIU=(PSIU-PSI)/DUT 
PHIV*(PHIV-PHI )/DUT 
PSIV=(PSIV-PSI)/DUT 

IF(IFACE.EQ.1 .AND.NOI .EQ.1)PSIU=PSIU-REC1/RB*2.D0*UB 
NMAX=2+I FACE 
LMAX=3- 1 FACE 
DO 5 N=1,NMAX 

C SET DISPLACEMENT AMPLITUDES 
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CALL ESET(ZB,N,E,DE) 

UR=0MDT*0I ( N )+UB/RB*T I ( N ) 

A(1,1)=(0.D0,0.D0) 

A(1,2)=(REC*(DV-VB/RB/HB*DRH)+PSIV+UR*REC1)*P1R 

A(1,3)=P1R*PSIU-TI(N)*REP*VB/RB 

A(2,1 )=(0.D0,0.D0) 

A(2,2)=DCMPLX(DRH/RB/HB,0.D0) 

A(2 < 3)=TI(N)/RB 

A(3,1)=TI(N)/RB/P1R 

A(3,2)=DCMPLX(DU*REC+PHIV,0.D0) 

A(3,3)=(PHIU+UR*REC1) 

C L=1 IS COMPLEMENTARY SOLUTION 
DO 7 1=1.3 

7 B(I,1)=(0.D0,0.D0) 

DO 6 L=2,LMAX 

TMP=(VB*(RB*DE(L)+IFACE*E(L))+RB*E(L)*(DV+UR))/RB/HB 

B( 1 , L)=P1R*PS I H*E ( L ) - REP*VB*TMP 

B(2,L)=TMP 

B(3,L)=DCMPLX<E(L)*PHIH,0.D0) 

6 CONTINUE 

C IF AXIAL (RADIAL) INERTIA IS INCLUDED DIVIDE U EQUATION BY COEFF OF DU/DS 
IF(NOI.NE.1)THEN 
NEQ=3 

DO 20 J=1 ,3 

20 A(3,J)=A(3,J)/REC/VB 
DO 21 L=1,LMAX 

21 B(3,L)=B(3,L)/REC/VB 
ELSE 

NEQ=2 

DO 23 1=1,2 
A( I ,3)=A( I ,3)/A(3,3) 

DO 24 L=1 ,LMAX 

24 B(I,L)=B(I,L)-A(I,3)*B(3,L) 

DO 23 J=1,2 

23 A(I,J)=A(I,J)-A(I,3)*A(3, J) 

END IF 

C REPLACE CB> WITH CB>- [AHY> 

DO 22 L=1 ,LMAX 
DO 22 1=1 ,NEQ 
DO 22 J=1 ,NEQ 

22 B(I ,L)=B(I ,L)-A(I , J)*Y(J,L,N) 

C REPLACE [AJ WITH [I]+DZ/2*(A] 

DO 8 1=1, NEQ 
DO 8 J=1 ,NEQ 
A(I , J)=DZ2*A(I, J) 

IF(I.EQ.J)A(I,J)=1.D0+A(I,J) 

8 CONTINUE 

C SOLVE EQUATIONS FOR ALL LMAX RIGHT HAND SIDE VECTORS IN ONE SHOT ■ 

CALL CMATIN(A,B, DETER, NEQ, LMAX, ID, 3, LABEL) 

1F(ID.NE.1)THEN 
IER=3 
RETURN 
END IF 

C CALCULATE NEW (Y> 

DO 9 L=1 ,LMAX 
DO 10 1=1, NEQ 

10 Y(I,L,N)=Y(I,L,N)+DZ*B(I,L) 

C UPDATE FORCE AND MOMENT INTEGRALS 

DFI=(Y(1,L,N)-DZ2*B(1,L))*RB*ABS(DZ) 

FOI(L,N)=FOI(L,N)+DFI 

IF(N.LT.3)EMI(L,N)=EMI(L,N)+ZB*DFI 

9 CONTINUE 

5 CONTINUE 

RETURN 

END 


SUBROUTINE ESET(Z,N,E,DE) 

C SETS DISPLACEMENT/TILT AMPLITUDE, E AND SLOPE DE 

C L=1 ,2,3 CORRESPONDS TO COMPLIMENTARY SOLUTION, TILT, RADIAL DISP RESP. WHEN N<3 

C L=1 ,2 CORRESPONDS TO COMP. SOL. AND AXIAL DISP RESP. WHEN N=3 
IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION E(3),DE(3),E0(3),DE0(3) 

DATA E0,DE0/0.D0,1 .00,1 .D0,0.D0,0.D0,0.D0/ 

DO 5 L=1,3 
DE(L)=DE0(L) 

5 E(L)=E0(L) 

IF(N.EQ.3)RETURN 
E(2)=Z 
DE(2)*1 .DO 
RETURN 


END 
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SUBROUT I NE KBCAKNOI , I FACE , I D I R , NREG, NRSUB , 

+RE,REC,P1R,0MT,0MDT1 ,DUT,ZET, 

+ 1 GROT , ALP , SBET , CBET , DELT , ENGP , ZETG ,UHG , VHG , 

+RG,ZTG,H,U,V,CK,CB,IER) 

C SETS UP BOUNDARY AND CONTINUITY CONDITIONS, SOLVES SECONDARY FLOW PROBLEM 
C AND CALCULATES STIFFNESS AND DAMPING COEFFICIENTS. 

C V (VI) 

C CALLED BY TSEAL 
C CALLS DELTP, ESET AND DSOLV 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

PARAMETER (NDZ=201,NDREG=21) 

C0MPLEX*16 Y(3,3,3),FOI(3,3),EMI(3,3) 

DIMENSION ALP(NDREG),SBET(NDREG),CBET(NDREG),DELT(NDREG), 
+IGROT(NDREG),UHG(NDZ,NDREG),VHG(NDZ,NDREG) 

DIMENSION NRSUB(NDREG),RG(NDZ,NDREG),ZET(NDREG),ZTG(NDZ,NDREG) 

DIMENSION H(NDZ,NDREG),U(NDZ,NDREG),V(NDZ,NDREG) 

DIMENSION ENGP(NDREG),ZETG(NDREG) 

DIMENSION E(3),DE(3),Y20(3),CK(4,4),CB(4,4) 

C INITIAL DISURBANCES IN INLET VELOCITY (V) FOR COMP AND PARTICULAR SOLUTIONS 
DATA Y20/1. DO, 0.00,0. DO/ 

RV20= . 5D0*REC*P1 R 
OMDT=OMDT1 

C AVOID INDETERMINACY FOR 0 FREQYENCY DISTURBANCE 

IF(ABS(0MDT).LT.1.D-4)0MDT=1.D-4 „. M1 

C SET UP STARTING CONDITIONS AND LOOPING PARAMETERS FOR SECONDARY FLOW SOLUTION 
KST=1 
KEN=NREG 

I F( IDIR.EQ.-1 JTHEN 
KST =NREG 
KEN=1 
END IF 

NMAX=2+IFACE 
LMAX=3-IFACE 
DO 30 K=KST,KEN, ID IR 
JST=1 

JEN=NRSUB(K) 

IF( IDIR. EQ .-1JTHEN 
JST=JEN+1 
JEN-2 
END IF 

HJMP=1 .DIO 

IF(K.NE.KST)HJMP=H( J1,K1)-H(JST,K) 

C GET DERIVATIVES FOR FLOW LOSS AT JUMP 
I FCNOI .NE.1)THEN 

CHI=DELTP(RE,V(JST,K),H(JST,K),HJMP,ZET(K>) „„, W „ 11T 

CHIH=(DELTP(RE,V(JST,K),H(JST,K)+OUT,HJMP,ZET(K))-CHI)/DUT 

CHIV=(DELTP(RE,V(JST,K)+DUT,H(JST,K),HJMP,ZET(K»-CHI)/DUT 

END IF 

DO 5 N=1,NMAX 

C SET DISPLACEMENT AMPLITUDES 

CALL ESET(ZTG(JST,K),N,E,DE) 

DO 5 L=1 ,LMAX 
IF(K.EQ.KST)THEN 

C SET UP INITIAL OR CONTINUITY CONDITIONS AT START OF EACH REGION 
C {Y> ARE DISTURBANCES IN PRESSURE, AXIAL VELOCITY AND TANGENTIAL VELOCITY 
Y(1 ,L,N)=(O.DO,O.DO) 

Y(2,L,N)=DCMPLX(Y20(L),0.D0) 

Y(3,L,N)=(0.D0,0.D0) 

FOI(L,N)=(O.DO,O.DO) 

EMI(L,N)=(0.D0,0.D0) 

ELSE 

Y(2,L,N)=V(JST,K)*HJMP/H(JST,K)/H(J1 ,K1)*E(L)+ 

+ H(J1 ,K1 )/H( JST,K)*Y(2,L,N) 

IF(NOI.NE.1)Y(1,L,N)=Y(1,L,N)+RV20*(-CHIH*E(L)+CHIV*Y(2,L,N)) 

5 CONTINUE 


K1 S K 

C STEP THROUGH REGION 

DO 30 J=J ST, JEN, IDIR 
J1=J+IDIR 

CALL DSOLV(NOI , I FACE ,RE,REC,P1R,OMT, OMDT , DUT , 

+ ALP(K),SBET(K),CBET(IO,DELT(K),ENGP(K),ZETG(IC), IGROT(K), 

+ RG( J.K),ZTG( J,K),H(J,K),U( J,K),V( J,K),UHG( J,IC),VHG(J,K), 

♦ RG(J1,K),ZTG(J1,K),H(J1,K),U(J1,K),V(J1,K), 

+ Y,FOI,EMI,IER) 

I F( IER.NE . 0 ) RETURN 
30 CONTINUE 

C COMBINE COMP AND PARTICULAR SOL. TO SATISFY P=0 AT DOWNSTREAM BOUNDARY 
DO 40 N=1,NMAX 
DO 40 L S 2,LMAX 

FOI(L,N)=-Y(1,L,N)/Y(1,1,N)*FOI<1,N)+FOI(L,N) 

40 EMI(L,N)*-Y(1,L,N)/Y(1,1,N)*EMI(1,N)+EMI(L,N) 

P1*4.D0*ATAN(1.D0) 

PI2=PI/2.D0 

C * EXTRACT " I? i FFN eIs AND DAMPING COEFFICIENTS FOR FACE SEAL 

C INITIALIZE STIFFNESS AND DAMPING MATRICES 
DO 41 1=1 ,4-IFACE 
DO 41 J=1 ,4-IFACE 
CK(I,J)=0.D0 

41 CB(I, J)*0.D0 

C AXIAL FORCE DUE TO AXIAL DISPLACEMENT 
CK(1,1)*2.D0*PI*DREAL(F0I(2,3)) 
CB(1,1)*2.DO*PI*DIMAG(FOI(2,3))/OMDT 
C MOMENTS DUE TO TILT 

CK(3,3)*PI2*DREAL(EMI(2,1)+EMI(2,2)) 

CK(2,2)=CK(3,3) 

CB(3,3)=PI2*DIMAG(EMI(2,1)-EMI(2,2))/0MDT 

CB(2,2)=CB(3,3) 

CK(2,3)=PI2*DIMAG(EMI(2,1)+EMI(2,2)) 

CK(3,2)=-CK(2,3) 

CB(2,3)=-PI2*DREAL(EMI(2,1)-EMI(2,2))/OMDT 

CB(3,2)*-CB(2,3) 

ELSE 

C EXTRACT STIFFNESS AND DAMPING COEFFICIENTS FOR CYLINDRICAL SEAL 
C MOMENTS DUE TO TILT 

CK<4,4)=PI2*DREAL<EMI(2,1)+EMI(2,2)> 

CK(3,3)*CK(4,4) 

CB(4,4)=PI2*DIMAG(EMI(2,1)-EMI (2,2))/OMDT 
CB(3,3)=CB<4,4) 

CK(3,4)=PI2*DIMAG(EMI(2,1)+EMI(2,2)) 

CK(4,3)=-CK(3,4) 

CB(3,4)=-PI2*DREAL(EMI(2,1)-EMI(2,2))/OMDT 

CB(4,3)=-CB(3,4) 

C MOMENTS DUE TO DISPLACEMENT 

CK(4,1 )=PI2*0REAL(EMI(3,1 )+EMI(3,2)) 

CK(3,2)*-CK(4,1) 

CB(4,1 )*PI2*DIMAG(EMI(3,1)-EMI(3,2))/0MDT 
CB(3,2)*-CB(4,1) 

CK(3,1)=PI2*DIMAG(EMI(3,1)+EMI(3,2» 

CK(4,2)-CK(3,1) 

CB(3,1)=-PI2*DREAL(EMI<3,1)-EMI(3,2))/OMDT 

CB(4,2)-CB(3,1) 

C FORCES DUE TO TILT 

CK(1,4)*PI2*DREAL(F0I(2,1)+F0I(2,2» 

CK(2,3)*-CK(1,4) 

CB(1,4)“PI2*DIMAG(FOI(2,1)-FOI(2,2))/OMDT 

CB(2,3)**CB(1,4> 

CK(2,4)*-PI2*DIMAG(FOI(2,1)+FOI(2,2)) 

CK(1 ,3) S CK(2,4) 

CB(2,4)=PI2*DREAL(FOI(2,1)-FOI(2,2))/OMDT 

CB(1,3)-CB(2,4) 
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C FORCES DUE TO DISPLACEMENT 

CK(1 ,1)=PI2*DREAL(FOI(3,1 )+F0I(3,2)) 
CK(2,2) S CK(1,1) 

CB(1,1)=PI2*DIMAG(FOI<3,1)-FOI(3,2))/OMDT 

CB(2,2)-CB(1,1) 

CK(2,1)— PI2*DIMAG(FOI(3,1)+FOI(3,2)) 

CK(1 ,2)— CK(2,1) 

CB(2,1)=PI2*DREAL(FOI(3,1)-FOI(3,2))/OMDT 
CB(1,2)— CB(2,1) 

END IF 

RETURN 

END 


to 

o 


SUBROUTINE CMATINCA, B, DETER, N1 , Ml , ID, N2, INDEX) 

C COMPEX MATRIX INVERTER 
C CALLED BY COLPC 

IMPLICIT C0MPLEX*16(A-H,0-Z) 

DOUBLE PRECISION AMAX 

DIMENSION A(N2,N2),B(N2,1), INDEX(N2,3) 

EQUIVALENCE <IR0W,JR0W), (ICOLU, JCOLU), (AMAX, T, SWAP) 

M=M1 

N=n 1 

10 DETER =(1 .D0,0.D0) 

DO 20 J=1 ,N 
20 INDEX(J,3) = 0 
DO 550 1-1, N 
AMAX-0.0D0 
DO 105 J=1 ,N 

IF(INDEX(J,3)-1) 60, 105, 60 
60 DO 100 K=1 ,N 

IF(INDEX(K,3)-1) 80, 100, 715 
80 IF ( AMAX -ABS (A(J,K))) 85, 100, 100 
85 IROW-J 
ICOLU-K 

AMAX -ABS (A(J,K)) 

100 CONTINUE 
105 CONTINUE 

I F(AMAX)1 10,715, 110 

110 INDEX(IC0LU,3) - INDEX(ICOLU,3) +1 
INDEX(I,1)-IROU 
INDEX(I,2)-IC0LU 

130 IF (I ROW* ICOLU) 140, 310, 140 
140 DETER— DETER 
DO 200 L*1 ,N 
SWAP=A(IROW,L) 

A(IROW,L)-A(ICOLU,L) 

200 A(ICOLU,L)=SWAP 
I F CM ) 310. 310, 210 
210 DO 250 L-1, M 
SWAP=B(IROW,L) 

B(IROW,L)=B(ICOLU,L) 

250 B( ICOLU, L)=SWAP 
310 PIVOT »A( ICOLU, ICOLU) 

IF(PIVOT.EQ.(O.DO,O.DO))GO TO 715 
DETER=DETER*PIVOT 
A(ICOLU,ICOLU)=(1 .DO.O.DO) 

DO 350 L=1,N 

350 AC I COLU , L )-A( I COLU , L )/P I VOT 
I F CM) 380, 380, 360 
360 DO 370 L=1,M 
370 B( ICOLU, L)=B( ICOLU, D/PIVOT 
380 DO 550 L1=1,N 

I FC LI - 1 COLU) 400, 550, 400 
400 T=A(L1, ICOLU) 

A(L1, ICOLU)*(O.DO,O.DO) 

IF(T.EQ.(O.DO,O.DO))GO TO 550 
430 DO 450 L=1,N 
450 A(L1,L)*A(L1,L)-A(ICOLU,L)*T 
IFCM) 550. 550, 460 
460 DO 500 L=1,M 
500 B( LI, L)=B(L1,L)-B( ICOLU, L)*T 
550 CONTINUE 
600 DO 710 1-1, N 

L=N+1 - I 

IF <INDEX(L,1)-INDEX(L,2)) 630, 710, 630 
630 JROU-INDEX(L,1) 

JCOLU=INDEX(L,2) 

DO 705 K-1.N 
SWAP=A(K, JROW) 
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A(K,JR0W)=A(K, JCOLU) 
A(K,JCOLU)=SWAP 
705 CONTINUE 
710 CONTINUE 
ID =1 

740 RETURN 
715 ID =2 

DETER=(0.D0,0.D0) 

RETURN 

END 


SUBROUTINE HATINV(A,B,DETER,N1 ,M1 , ID,N2, INDEX) 

C REAL MATRIX INVERSION ROUTINE 
C CALLED BY HOME COLP 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION A(N2,N2),B(N2,1 ), INDEX(N2,3) 
EQUIVALENCE <IROW,JROW), (ICOLU, JCOLU), (AMAX, T 
M=M1 
N=N1 

10 DETER * 1.D0 
DO 20 J=1,N 
20 INDEX(J,3) = 0 
DO 550 1=1 ,N 
AMAX=0.D0 
DO 105 J=1 ,N 

IF(INDEX(J,3)-1) 60, 105, 60 
60 DO 100 K=1,N 

IF(INDEX(K,3)-1) 80, 100, 715 
80 IF ( AMAX -ABS (A(J,K))) 85, 100, 100 
85 IROW=J 
ICOLU-K 

AMAX = ABS (A(J,K)) 

100 CONTINUE 
105 CONTINUE 

IF(AMAX)110,715,110 

110 INDEX(ICOLU,3) = INDEX(IC0LU,3) +1 
INDEX(I,1)=IROW 
INDEX(I,2)=ICOLU 

130 IF (I ROW- ICOLU) 140, 310, 140 
140 DETER=-DETER 
DO 200 L=1 ,N 
SWAP=A( I ROW, L ) 

A(IROW,L)=A(ICOLU,L) 

200 A(ICOLU,L)=SWAP 
IF(M) 310. 310, 210 
210 DO 250 L=1, M 
SWAP=B(IROU,L) 

B(IROW,L)=B(ICOLU,L) 

250 B(ICOLU,L)=SWAP 
310 PIVOT =A( ICOLU, ICOLU) 

I F(PIVOT.EQ.O.DO)GO TO 715 
DETER=DETER*PIVOT 
A( ICOLU, ICOLU)=1 .DO 
DO 350 L=1 ,N 

350 A( ICOLU, L)=A( ICOLU, D/PIVOT 
IF(M) 360, 380, 360 
360 DO 370 L=1 ,M 
370 B( ICOLU, L)=B( ICOLU, D/PIVOT 
380 DO 550 L1=1,N 

IF(LI-ICOLU) 400, 550, 400 
400 T=A(L1, ICOLU) 

ACL1 , ICOLU)=O.DO 
IF(T)430,550,430 
430 DO 450 L=1,N 

450 A(L1,L)=A(L1 ,L)-A( ICOLU, L)*T 
IF(M) 550, 550, 460 
460 DO 500 L=1 ,M 
500 B(L1,L)=B(L1 ,L)-B( ICOLU, L)*T 
550 CONTINUE 
600 DO 710 1=1, N 
L=N+1 - 1 

IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630 
630 JROU=INDEX(L,1) 

J COLU= I NDEX( L , 2 ) 

DO 705 K=1 ,N 
SWAP=A(K, JROW) 

A(K,JROW)=A(K, JCOLU) 
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A(K, JCOLU) s SWAP 
705 CONTINUE 
710 CONTINUE 
ID =1 

740 RETURN 
715 ID =2 

DETER=0.D0 

RETURN 

END 
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